library(sf)
library(mapview)
library(dplyr)
library(tidyr)
library(lubridate)
library(evd)
library(parallel)
library(FNN)
library(stringr)
pathshp <- "~/comp9991/data/geo/DOPGrid.shp"
map <- st_read(pathshp, quiet = TRUE)
print(map)
Simple feature collection with 17450 features and 4 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 274189.9 ymin: 2901792 xmax: 532797 ymax: 3324561
Projected CRS: NAD83 / UTM zone 17N
First 10 features:
OBJECTID PIXEL Shape_Leng Shape_Area geometry
1 14202 131641 7697.982 3699323 POLYGON ((411365.9 3234204,...
2 14203 131642 7697.949 3699292 POLYGON ((413224.1 3234190,...
3 14204 131643 7697.920 3699264 POLYGON ((415082.4 3234176,...
4 14205 131644 7697.892 3699237 POLYGON ((416940.6 3234162,...
5 14206 131645 7697.862 3699208 POLYGON ((418798.9 3234149,...
6 14207 131646 7697.834 3699181 POLYGON ((420657.1 3234136,...
7 14208 131647 7697.806 3699155 POLYGON ((422515.3 3234123,...
8 14209 131648 7697.780 3699130 POLYGON ((424373.5 3234111,...
9 14210 131649 7697.756 3699106 POLYGON ((426231.8 3234098,...
10 14211 131650 7697.728 3699080 POLYGON ((428090 3234087, 4...
mapview(map)
data <- read.csv("~/comp9991/data/PixelRain15min_1995.csv",header = TRUE)
change_time_resolution <- function(data, res = "default"){
# data: dataframe
data <- data %>%
mutate(DateTime = dmy_hms(DateTime),
Date = as.Date(DateTime),
Hour = floor_date(DateTime, "hour"))
if (res == "default"){
return(data)
}
else if (res == "hourly"){
data_hourly <- data %>%
group_by(Hour) %>%
summarise(across(starts_with("P"), \(x) sum(x, na.rm = TRUE)))
return(data_hourly)
}
else if (res == "daily"){
data_daily <- data %>%
group_by(Date) %>%
summarise(across(starts_with("P"), \(x) sum(x, na.rm = TRUE)))
return(data_daily)
}
}
hourly
data_hourly <- change_time_resolution(data,res = "hourly")
class(data_hourly)
[1] "tbl_df" "tbl" "data.frame"
all_pixel_IDs <- names(data_hourly)[-1]
all_pixel_IDs <- as.numeric(stringi::stri_extract_first(all_pixel_IDs, regex = "[0-9]+"))
sub_map <- map %>%
filter(PIXEL %in% all_pixel_IDs)
mapview(sub_map)
cut a square: 14*14
cutid1 <- seq(11951+2, 11951 + 13+2)
cutid2 <- seq(11874+2, 11874 + 13+2)
cutid3 <- seq(11796+2, 11796 + 13+2)
cutid4 <- seq(11716+2, 11716 + 13+2)
cutid5 <- seq(11635+2, 11635 + 13+2)
cutid6 <- seq(11553+2, 11553 + 13+2)
cutid7 <- seq(11470+2, 11470 + 13+2)
cutid8 <- seq(11386+2, 11386 + 13+2)
cutid9 <- seq(11301+2, 11301 + 13+2)
cutid10 <- seq(11215+2, 11215 + 13+2)
cutid11 <- seq(11129+2, 11129 + 13+2)
cutid12 <- seq(11042+2, 11042 + 13+2)
cutid13 <- seq(10954+2, 10954 + 13+2)
cutid14 <- seq(10866+2, 10866 + 13+2)
cut_id <- c(
cutid1,
cutid2,
cutid3,
cutid4,
cutid5,
cutid6,
cutid7,
cutid8,
cutid9,
cutid10,
cutid11,
cutid12,
cutid13,
cutid14
)
square_map <- sub_map %>%
filter(OBJECTID %in% cut_id)
mapview(square_map)
data_hourly <- data_hourly[rowSums(data_hourly[, -1]) > 0, ]
sub_pixel <- square_map$PIXEL
idx <- all_pixel_IDs %in% sub_pixel
idx <- c(TRUE, idx)
data_hourly_sub <- data_hourly[, idx]
val <- data_hourly_sub[20, 2:ncol(data_hourly_sub)]
val_long <- val %>%
pivot_longer(
cols = starts_with("P"),
names_to = "PIXEL_COL_NAME",
values_to = "value"
) %>%
mutate(
PIXEL = as.integer(str_remove(PIXEL_COL_NAME, pattern = "P"))
) %>%
select(PIXEL, value)
square_map_value <- square_map %>%
left_join(val_long, by = "PIXEL")
mapview(square_map_value, zcol="value")
filter extremes
rain_hourly_mat <- as.matrix(data_hourly_sub[, -1])
rain_hourly_uniform <- apply(rain_hourly_mat, 2, function(col) {
ind <- col > 0
u <- numeric(length(col))
u[ind] <- rank(col[ind]) / (sum(ind) + 1)
u
})
eps <- 1e-10
rain_hourly_uniform[rain_hourly_uniform <= 0] <- eps
rain_hourly_uniform[rain_hourly_uniform >= 1] <- 1 - eps
rain_hourly_pareto <- qgpd(rain_hourly_uniform, xi = 1, beta = 1, mu = 1)
L1_risk_hourly <- rowMeans(rain_hourly_pareto, na.rm = TRUE)
u_L1_hourly <- quantile(L1_risk_hourly, 0.90, na.rm = TRUE)
extreme_L1_hourly <- rain_hourly_pareto[L1_risk_hourly > u_L1_hourly, ]
idx of timestep
idx_ext_hourly <- which(L1_risk_hourly > u_L1_hourly & !is.na(L1_risk_hourly))
times_hourly <- data_hourly_sub[[1]][idx_ext_hourly]
square_map_m <- st_transform(square_map, crs = 3577)
centroids_m <- st_centroid(square_map_m)
Warning: st_centroid assumes attributes are constant over geometries
coords_m <- st_coordinates(centroids_m)
nn <- get.knn(coords_m, k = 2)$nn.dist[, 2]
cell_len_m <- median(nn)
coord <- coords_m / cell_len_m
build_Gamma <- function(coord, alpha, theta){
D <- as.matrix(dist(coord))
G <- (D / theta)^alpha
diag(G) <- 0 # \gamma(0) = 0
return(G)
}
prep_br <- function(Gamma) {
d <- nrow(Gamma)
out <- vector("list", d)
for (j in 1:d) {
idx <- setdiff(seq_len(d), j)
Sigj <- ( outer(Gamma[idx, j], rep(1, d-1)) +
outer(rep(1, d-1), Gamma[idx, j]) -
Gamma[idx, idx] )
Sigj <- Sigj + .Machine$double.eps * diag(d-1)
Lj <- chol(Sigj)
out[[j]] <- list(idx = idx, L = Lj, Gamma_col = Gamma[, j])
}
out
}
intensity_logskew <- function(x,par,alpha.para=TRUE,log=TRUE){
sigma = par[[1]]
if(!is.matrix(x)){x <- matrix(x,nrow=1)}
n = ncol(x)
if(n==1) return(1/(x^2))
omega2 = diag(sigma)
chol.sigma = chol(sigma)
inv.sigma = chol2inv(chol.sigma)
logdet.sigma = sum(log(diag(chol.sigma)))*2
if(alpha.para){
alpha = par[[2]]
delta = c(sigma %*% alpha)/sqrt(c(1+alpha %*% sigma %*% alpha))
}else{
delta = par[[2]]
alpha = c(1 - delta %*% inv.sigma %*% delta)^(-1/2) * c(inv.sigma %*% delta)
}
a = log(2) + pnorm(delta,log.p=TRUE)
q = rowSums(inv.sigma)
sum.q = sum(q);sum.alpha = sum(alpha)
q.mat = matrix(q,n,n,byrow=TRUE)
x.log = log(x)
x.circ = x.log + matrix(a,nrow=nrow(x),ncol=n,byrow=TRUE)
beta = (1+sum.alpha^2/sum.q)^(-0.5)
tau.tilde = apply(x.circ,1,function(x.i) beta * sum((alpha - sum.alpha*q/sum.q) * (x.i + omega2/2))+ beta*sum.alpha/sum.q)
A = inv.sigma - q %*% t(q)/sum.q
val = -(n-3)/2 * log(2) - (n-1)/2*log(pi)-1/2*logdet.sigma - 1/2*log(sum.q) - 1/2 * (sum(q*omega2)-1)/sum.q - 1/8*c(omega2 %*% A %*% omega2)
val = val - rowSums(x.log) - 1/2 * apply(x.circ,1,function(x.i) c(x.i %*% A %*% x.i) + sum(x.i * (2*q/sum.q + c(A %*% omega2)))) + pnorm(tau.tilde,log.p=TRUE)
if(log)
return(val)
else
return(exp(val))
}
dlog_beta <- function(x, a, b) dbeta(x, a, b, log = TRUE) # mixture weight prior
dlog_gamma <- function(x, a, b) dgamma(x, shape = a, rate = b, log = TRUE)
component_loglik <- function(X, coord, alpha, theta, anchor = 1){
Gamma <- build_Gamma(coord = coord, alpha = alpha, theta = theta)
prep <- prep_br(Gamma)
idx <- prep[[anchor]]$idx
L <- prep[[anchor]]$L
Sigma <- t(L) %*% L
Xsub <- as.matrix(X[, idx, drop = FALSE])
intensity_logskew(
Xsub, par = list(Sigma, delta = rep(0, ncol(Xsub))),
alpha.para = FALSE, log = TRUE
)
}
mix_loglik <- function(X, coord, alpha, th1, th2, anchor = 1) {
if (!is.finite(th1) || th1 <= 0) return(list(l1 = rep(-Inf, nrow(X)), l2 = rep(-Inf, nrow(X))))
if (!is.finite(th2) || th2 <= 0) return(list(l1 = rep(-Inf, nrow(X)), l2 = rep(-Inf, nrow(X))))
l1 <- component_loglik(X, coord, alpha = alpha, theta = th1, anchor = anchor)
l2 <- component_loglik(X, coord, alpha = alpha, theta = th2, anchor = anchor)
list(l1 = l1, l2 = l2)
}
tau_from_ll <- function(l1, l2, lambda) {
a <- log(lambda) + l1
b <- log(1 - lambda) + l2
m <- pmax(a, b)
exp(a - m) / (exp(a - m) + exp(b - m))
}
mcmc_br <- function(
data, coord, alpha,
n_iter = 4000,
burn = 1000,
thin = 2,
# priors
a_lambda = 1, b_lambda = 1, # for mixture weight lambda
a_th2 = 1, b_th2 = 1, # for theta2
a_del = 1, b_del = 1, # for delta
# proposal sds on log-scale
prop_sd_log_th2 = 0.10,
prop_sd_log_del = 0.10,
# initials
init_lambda = 0.5, # mixture weight
init_th2 = 15, # theta2
init_del = 20, # delta
anchor = 1,
verbose_every = 100
){
stopifnot(init_lambda > 0, init_lambda < 1, init_th2 > 0, init_del > 0)
keep_idx <- seq.int(burn + 1, n_iter, by = thin)
n_keep <- length(keep_idx)
draws <- data.frame(
iter = keep_idx,
lambda = NA_real_, # mixture weight
theta1 = NA_real_,
theta2 = NA_real_
)
z_draws <- matrix(NA_integer_, nrow = n_keep, ncol = nrow(data))
tau_bar <- rep(0, nrow(data))
# init
lambda <- init_lambda
th2 <- init_th2
del <- init_del
th1 <- th2 + del
ll <- mix_loglik(X = data, coord = coord, alpha = alpha, th1 = th1, th2 = th2, anchor = anchor)
tau <- tau_from_ll(ll$l1, ll$l2, lambda = lambda)
z <- rbinom(n = length(tau), size = 1, prob = tau)
acc_block <- 0L; n_prop <- 0L
# helper: complete-data log-posterior in η-space (η2=log θ2, ηd=log δ)
logpost_eta <- function(eta2, etad, z_vec) {
th2p <- exp(eta2)
delp <- exp(etad)
th1p <- th2p + delp
llp <- mix_loglik(X = data, coord = coord, alpha = alpha,
th1 = th1p, th2 = th2p, anchor = anchor)
lp_ll <- sum(z_vec * llp$l1 + (1 - z_vec) * llp$l2)
lp_prior <- dlog_gamma(th2p, a_th2, b_th2) + dlog_gamma(delp, a_del, b_del)
lp_jac <- eta2 + etad
list(value = lp_ll + lp_prior + lp_jac,
th1 = th1p, th2 = th2p, del = delp, l1_vec = llp$l1, l2_vec = llp$l2)
}
keep_cursor <- 1L
for (it in seq_len(n_iter)) {
## (1) blocked MH for (θ2, δ) on log-scale
n_prop <- n_prop + 1L
eta2 <- log(th2); etad <- log(del)
eta2p <- rnorm(1, mean = eta2, sd = prop_sd_log_th2)
etadp <- rnorm(1, mean = etad, sd = prop_sd_log_del)
# current value via same function
curr <- logpost_eta(eta2, etad, z)
prop <- logpost_eta(eta2p, etadp, z)
if (log(runif(1)) < (prop$value - curr$value)) {
th2 <- prop$th2; del <- prop$del; th1 <- prop$th1
ll$l1 <- prop$l1_vec; ll$l2 <- prop$l2_vec
acc_block <- acc_block + 1L
}
# else keep (th1,th2,del,ll) as-is
## (2) Gibbs for lambda | z (Beta-conjugate)
n1 <- sum(z); n0 <- length(z) - n1
lambda <- rbeta(1, a_lambda + n1, b_lambda + n0)
## (3) Gibbs for z | x, θ, λ
tau <- tau_from_ll(ll$l1, ll$l2, lambda)
z <- rbinom(length(tau), 1, tau)
## store thinned draws
if (it %in% keep_idx) {
tau_bar <- tau_bar + tau / n_keep
draws$lambda[keep_cursor] <- lambda
draws$theta1[keep_cursor] <- th1
draws$theta2[keep_cursor] <- th2
z_draws[keep_cursor, ] <- z
keep_cursor <- keep_cursor + 1L
}
if (verbose_every > 0 && it %% verbose_every == 0) {
cat(sprintf("iter %d | lambda=%.3f θ1=%.3f θ2=%.3f | acc(θ-block)=%.3f\n",
it, lambda, th1, th2, acc_block / n_prop))
}
}
list(
draws = draws,
z_draws = z_draws,
accept_rate_theta_block = acc_block / n_prop,
tau_bar = tau_bar,
hard_labels = ifelse(tau_bar > 0.5, 1, 2),
post_mean = c(lambda = mean(draws$lambda),
theta1 = mean(draws$theta1),
theta2 = mean(draws$theta2)),
post_median = c(lambda = median(draws$lambda),
theta1 = median(draws$theta1),
theta2 = median(draws$theta2))
)
}
fit the data
alpha <- 1.5
fit1 <- mcmc_br(
data = extreme_L1_hourly, coord = coord, alpha = alpha,
n_iter = 10000, burn = 2000, thin = 2,
# priors
a_lambda = 1, b_lambda = 1, # prior for mixture weight lambda
a_th2 = 1, b_th2 = 1, # prior for theta2
a_del = 1, b_del = 1, # prior for delta
# initials
init_th2 = 2,
init_del = 4,
# proposals (log-scale)
prop_sd_log_th2 = 0.01,
prop_sd_log_del = 0.01,
anchor = 1,
verbose_every = 50
)
iter 50 | lambda=0.429 θ1=6.866 θ2=2.087 | acc(θ-block)=0.440
iter 100 | lambda=0.419 θ1=6.927 θ2=2.082 | acc(θ-block)=0.330
iter 150 | lambda=0.421 θ1=6.854 θ2=2.088 | acc(θ-block)=0.287
iter 200 | lambda=0.402 θ1=6.889 θ2=2.100 | acc(θ-block)=0.265
iter 250 | lambda=0.379 θ1=6.896 θ2=2.104 | acc(θ-block)=0.276
iter 300 | lambda=0.393 θ1=6.877 θ2=2.098 | acc(θ-block)=0.297
iter 350 | lambda=0.390 θ1=6.863 θ2=2.081 | acc(θ-block)=0.300
iter 400 | lambda=0.422 θ1=6.822 θ2=2.092 | acc(θ-block)=0.285
iter 450 | lambda=0.363 θ1=6.900 θ2=2.092 | acc(θ-block)=0.280
iter 500 | lambda=0.417 θ1=6.800 θ2=2.094 | acc(θ-block)=0.274
iter 550 | lambda=0.428 θ1=6.818 θ2=2.084 | acc(θ-block)=0.291
iter 600 | lambda=0.376 θ1=6.847 θ2=2.092 | acc(θ-block)=0.288
iter 650 | lambda=0.447 θ1=6.826 θ2=2.092 | acc(θ-block)=0.275
iter 700 | lambda=0.397 θ1=6.877 θ2=2.083 | acc(θ-block)=0.283
iter 750 | lambda=0.388 θ1=6.921 θ2=2.087 | acc(θ-block)=0.285
iter 800 | lambda=0.378 θ1=6.824 θ2=2.088 | acc(θ-block)=0.286
iter 850 | lambda=0.411 θ1=6.806 θ2=2.087 | acc(θ-block)=0.287
iter 900 | lambda=0.426 θ1=6.891 θ2=2.106 | acc(θ-block)=0.289
iter 950 | lambda=0.444 θ1=6.879 θ2=2.102 | acc(θ-block)=0.288
iter 1000 | lambda=0.428 θ1=6.795 θ2=2.092 | acc(θ-block)=0.288
iter 1050 | lambda=0.383 θ1=6.841 θ2=2.084 | acc(θ-block)=0.293
iter 1100 | lambda=0.426 θ1=6.824 θ2=2.083 | acc(θ-block)=0.292
iter 1150 | lambda=0.382 θ1=6.842 θ2=2.098 | acc(θ-block)=0.295
iter 1200 | lambda=0.419 θ1=6.878 θ2=2.094 | acc(θ-block)=0.293
iter 1250 | lambda=0.400 θ1=6.866 θ2=2.111 | acc(θ-block)=0.296
iter 1300 | lambda=0.393 θ1=6.793 θ2=2.089 | acc(θ-block)=0.293
iter 1350 | lambda=0.408 θ1=6.819 θ2=2.099 | acc(θ-block)=0.293
iter 1400 | lambda=0.413 θ1=6.850 θ2=2.090 | acc(θ-block)=0.286
iter 1450 | lambda=0.425 θ1=6.883 θ2=2.094 | acc(θ-block)=0.291
iter 1500 | lambda=0.418 θ1=6.905 θ2=2.086 | acc(θ-block)=0.293
iter 1550 | lambda=0.396 θ1=6.833 θ2=2.091 | acc(θ-block)=0.291
iter 1600 | lambda=0.415 θ1=6.790 θ2=2.091 | acc(θ-block)=0.292
iter 1650 | lambda=0.421 θ1=6.896 θ2=2.081 | acc(θ-block)=0.293
iter 1700 | lambda=0.404 θ1=6.905 θ2=2.093 | acc(θ-block)=0.291
iter 1750 | lambda=0.407 θ1=6.872 θ2=2.084 | acc(θ-block)=0.295
iter 1800 | lambda=0.429 θ1=6.812 θ2=2.103 | acc(θ-block)=0.297
iter 1850 | lambda=0.405 θ1=6.776 θ2=2.095 | acc(θ-block)=0.299
iter 1900 | lambda=0.431 θ1=6.843 θ2=2.088 | acc(θ-block)=0.300
iter 1950 | lambda=0.434 θ1=6.839 θ2=2.095 | acc(θ-block)=0.298
iter 2000 | lambda=0.416 θ1=6.843 θ2=2.085 | acc(θ-block)=0.296
iter 2050 | lambda=0.422 θ1=6.875 θ2=2.093 | acc(θ-block)=0.298
iter 2100 | lambda=0.401 θ1=6.827 θ2=2.089 | acc(θ-block)=0.298
iter 2150 | lambda=0.412 θ1=6.935 θ2=2.082 | acc(θ-block)=0.297
iter 2200 | lambda=0.426 θ1=6.834 θ2=2.096 | acc(θ-block)=0.298
iter 2250 | lambda=0.370 θ1=6.881 θ2=2.090 | acc(θ-block)=0.299
iter 2300 | lambda=0.465 θ1=6.816 θ2=2.096 | acc(θ-block)=0.297
iter 2350 | lambda=0.427 θ1=6.869 θ2=2.088 | acc(θ-block)=0.297
iter 2400 | lambda=0.450 θ1=6.830 θ2=2.095 | acc(θ-block)=0.296
iter 2450 | lambda=0.444 θ1=6.794 θ2=2.086 | acc(θ-block)=0.294
iter 2500 | lambda=0.412 θ1=6.773 θ2=2.082 | acc(θ-block)=0.294
iter 2550 | lambda=0.424 θ1=6.824 θ2=2.107 | acc(θ-block)=0.294
iter 2600 | lambda=0.387 θ1=6.920 θ2=2.109 | acc(θ-block)=0.293
iter 2650 | lambda=0.413 θ1=6.863 θ2=2.099 | acc(θ-block)=0.295
iter 2700 | lambda=0.444 θ1=6.824 θ2=2.091 | acc(θ-block)=0.295
iter 2750 | lambda=0.397 θ1=6.922 θ2=2.106 | acc(θ-block)=0.293
iter 2800 | lambda=0.375 θ1=6.832 θ2=2.100 | acc(θ-block)=0.297
iter 2850 | lambda=0.409 θ1=6.829 θ2=2.091 | acc(θ-block)=0.296
iter 2900 | lambda=0.438 θ1=6.822 θ2=2.081 | acc(θ-block)=0.298
iter 2950 | lambda=0.435 θ1=6.937 θ2=2.100 | acc(θ-block)=0.298
iter 3000 | lambda=0.472 θ1=6.796 θ2=2.085 | acc(θ-block)=0.299
iter 3050 | lambda=0.402 θ1=6.873 θ2=2.094 | acc(θ-block)=0.299
iter 3100 | lambda=0.437 θ1=6.870 θ2=2.094 | acc(θ-block)=0.298
iter 3150 | lambda=0.445 θ1=6.828 θ2=2.093 | acc(θ-block)=0.299
iter 3200 | lambda=0.420 θ1=6.910 θ2=2.088 | acc(θ-block)=0.299
iter 3250 | lambda=0.431 θ1=6.842 θ2=2.097 | acc(θ-block)=0.298
iter 3300 | lambda=0.411 θ1=6.885 θ2=2.088 | acc(θ-block)=0.299
iter 3350 | lambda=0.394 θ1=6.894 θ2=2.098 | acc(θ-block)=0.299
iter 3400 | lambda=0.395 θ1=6.814 θ2=2.107 | acc(θ-block)=0.299
iter 3450 | lambda=0.380 θ1=6.784 θ2=2.092 | acc(θ-block)=0.300
iter 3500 | lambda=0.484 θ1=6.853 θ2=2.087 | acc(θ-block)=0.300
iter 3550 | lambda=0.384 θ1=6.856 θ2=2.095 | acc(θ-block)=0.300
iter 3600 | lambda=0.427 θ1=6.850 θ2=2.084 | acc(θ-block)=0.301
iter 3650 | lambda=0.404 θ1=6.902 θ2=2.097 | acc(θ-block)=0.300
iter 3700 | lambda=0.423 θ1=6.773 θ2=2.089 | acc(θ-block)=0.299
iter 3750 | lambda=0.459 θ1=6.865 θ2=2.110 | acc(θ-block)=0.300
iter 3800 | lambda=0.406 θ1=6.808 θ2=2.088 | acc(θ-block)=0.301
iter 3850 | lambda=0.414 θ1=6.800 θ2=2.097 | acc(θ-block)=0.300
iter 3900 | lambda=0.421 θ1=6.834 θ2=2.097 | acc(θ-block)=0.300
iter 3950 | lambda=0.432 θ1=6.845 θ2=2.104 | acc(θ-block)=0.300
iter 4000 | lambda=0.431 θ1=6.840 θ2=2.095 | acc(θ-block)=0.299
iter 4050 | lambda=0.412 θ1=6.800 θ2=2.090 | acc(θ-block)=0.300
iter 4100 | lambda=0.409 θ1=6.854 θ2=2.096 | acc(θ-block)=0.299
iter 4150 | lambda=0.392 θ1=6.789 θ2=2.090 | acc(θ-block)=0.300
iter 4200 | lambda=0.414 θ1=6.844 θ2=2.097 | acc(θ-block)=0.300
iter 4250 | lambda=0.434 θ1=6.853 θ2=2.077 | acc(θ-block)=0.300
iter 4300 | lambda=0.387 θ1=6.849 θ2=2.095 | acc(θ-block)=0.299
iter 4350 | lambda=0.414 θ1=6.901 θ2=2.097 | acc(θ-block)=0.297
iter 4400 | lambda=0.409 θ1=6.923 θ2=2.086 | acc(θ-block)=0.297
iter 4450 | lambda=0.413 θ1=6.781 θ2=2.091 | acc(θ-block)=0.297
iter 4500 | lambda=0.420 θ1=6.883 θ2=2.095 | acc(θ-block)=0.298
iter 4550 | lambda=0.437 θ1=6.876 θ2=2.091 | acc(θ-block)=0.299
iter 4600 | lambda=0.375 θ1=6.823 θ2=2.099 | acc(θ-block)=0.298
iter 4650 | lambda=0.438 θ1=6.866 θ2=2.091 | acc(θ-block)=0.298
iter 4700 | lambda=0.420 θ1=6.817 θ2=2.094 | acc(θ-block)=0.299
iter 4750 | lambda=0.412 θ1=6.921 θ2=2.094 | acc(θ-block)=0.299
iter 4800 | lambda=0.391 θ1=6.841 θ2=2.093 | acc(θ-block)=0.299
iter 4850 | lambda=0.402 θ1=6.868 θ2=2.100 | acc(θ-block)=0.299
iter 4900 | lambda=0.390 θ1=6.861 θ2=2.100 | acc(θ-block)=0.300
iter 4950 | lambda=0.407 θ1=6.825 θ2=2.090 | acc(θ-block)=0.300
iter 5000 | lambda=0.395 θ1=6.845 θ2=2.096 | acc(θ-block)=0.299
iter 5050 | lambda=0.407 θ1=6.833 θ2=2.090 | acc(θ-block)=0.298
iter 5100 | lambda=0.396 θ1=6.895 θ2=2.091 | acc(θ-block)=0.298
iter 5150 | lambda=0.368 θ1=6.906 θ2=2.109 | acc(θ-block)=0.297
iter 5200 | lambda=0.413 θ1=6.820 θ2=2.086 | acc(θ-block)=0.299
iter 5250 | lambda=0.426 θ1=6.832 θ2=2.099 | acc(θ-block)=0.298
iter 5300 | lambda=0.457 θ1=7.003 θ2=2.097 | acc(θ-block)=0.298
iter 5350 | lambda=0.407 θ1=6.853 θ2=2.087 | acc(θ-block)=0.298
iter 5400 | lambda=0.396 θ1=6.839 θ2=2.095 | acc(θ-block)=0.298
iter 5450 | lambda=0.373 θ1=6.909 θ2=2.092 | acc(θ-block)=0.298
iter 5500 | lambda=0.405 θ1=6.892 θ2=2.080 | acc(θ-block)=0.298
iter 5550 | lambda=0.404 θ1=6.831 θ2=2.089 | acc(θ-block)=0.298
iter 5600 | lambda=0.442 θ1=6.887 θ2=2.088 | acc(θ-block)=0.297
iter 5650 | lambda=0.386 θ1=6.843 θ2=2.078 | acc(θ-block)=0.297
iter 5700 | lambda=0.393 θ1=6.896 θ2=2.096 | acc(θ-block)=0.297
iter 5750 | lambda=0.396 θ1=6.903 θ2=2.101 | acc(θ-block)=0.297
iter 5800 | lambda=0.429 θ1=6.786 θ2=2.094 | acc(θ-block)=0.297
iter 5850 | lambda=0.437 θ1=6.881 θ2=2.083 | acc(θ-block)=0.297
iter 5900 | lambda=0.443 θ1=6.841 θ2=2.098 | acc(θ-block)=0.297
iter 5950 | lambda=0.409 θ1=6.806 θ2=2.096 | acc(θ-block)=0.297
iter 6000 | lambda=0.426 θ1=6.859 θ2=2.097 | acc(θ-block)=0.297
iter 6050 | lambda=0.379 θ1=6.868 θ2=2.092 | acc(θ-block)=0.298
iter 6100 | lambda=0.402 θ1=6.841 θ2=2.103 | acc(θ-block)=0.298
iter 6150 | lambda=0.427 θ1=6.837 θ2=2.093 | acc(θ-block)=0.299
iter 6200 | lambda=0.413 θ1=6.866 θ2=2.069 | acc(θ-block)=0.300
iter 6250 | lambda=0.425 θ1=6.815 θ2=2.083 | acc(θ-block)=0.300
iter 6300 | lambda=0.459 θ1=6.787 θ2=2.097 | acc(θ-block)=0.300
iter 6350 | lambda=0.415 θ1=6.884 θ2=2.107 | acc(θ-block)=0.300
iter 6400 | lambda=0.355 θ1=6.890 θ2=2.096 | acc(θ-block)=0.300
iter 6450 | lambda=0.413 θ1=6.852 θ2=2.079 | acc(θ-block)=0.300
iter 6500 | lambda=0.429 θ1=6.851 θ2=2.105 | acc(θ-block)=0.301
iter 6550 | lambda=0.449 θ1=6.868 θ2=2.075 | acc(θ-block)=0.300
iter 6600 | lambda=0.399 θ1=6.870 θ2=2.098 | acc(θ-block)=0.300
iter 6650 | lambda=0.432 θ1=6.878 θ2=2.107 | acc(θ-block)=0.300
iter 6700 | lambda=0.423 θ1=6.850 θ2=2.099 | acc(θ-block)=0.300
iter 6750 | lambda=0.437 θ1=6.833 θ2=2.099 | acc(θ-block)=0.299
iter 6800 | lambda=0.386 θ1=6.915 θ2=2.101 | acc(θ-block)=0.300
iter 6850 | lambda=0.443 θ1=6.844 θ2=2.093 | acc(θ-block)=0.300
iter 6900 | lambda=0.398 θ1=6.832 θ2=2.103 | acc(θ-block)=0.299
iter 6950 | lambda=0.385 θ1=6.863 θ2=2.100 | acc(θ-block)=0.299
iter 7000 | lambda=0.353 θ1=6.860 θ2=2.104 | acc(θ-block)=0.299
iter 7050 | lambda=0.455 θ1=6.841 θ2=2.089 | acc(θ-block)=0.299
iter 7100 | lambda=0.396 θ1=6.833 θ2=2.095 | acc(θ-block)=0.298
iter 7150 | lambda=0.457 θ1=6.860 θ2=2.080 | acc(θ-block)=0.299
iter 7200 | lambda=0.412 θ1=6.890 θ2=2.101 | acc(θ-block)=0.299
iter 7250 | lambda=0.403 θ1=6.835 θ2=2.089 | acc(θ-block)=0.299
iter 7300 | lambda=0.458 θ1=6.819 θ2=2.083 | acc(θ-block)=0.299
iter 7350 | lambda=0.385 θ1=6.919 θ2=2.083 | acc(θ-block)=0.300
iter 7400 | lambda=0.426 θ1=6.829 θ2=2.109 | acc(θ-block)=0.299
iter 7450 | lambda=0.454 θ1=6.874 θ2=2.100 | acc(θ-block)=0.299
iter 7500 | lambda=0.460 θ1=6.803 θ2=2.097 | acc(θ-block)=0.299
iter 7550 | lambda=0.413 θ1=6.808 θ2=2.104 | acc(θ-block)=0.298
iter 7600 | lambda=0.453 θ1=6.820 θ2=2.092 | acc(θ-block)=0.298
iter 7650 | lambda=0.386 θ1=6.839 θ2=2.098 | acc(θ-block)=0.298
iter 7700 | lambda=0.418 θ1=6.844 θ2=2.098 | acc(θ-block)=0.299
iter 7750 | lambda=0.376 θ1=6.866 θ2=2.082 | acc(θ-block)=0.300
iter 7800 | lambda=0.437 θ1=6.842 θ2=2.083 | acc(θ-block)=0.300
iter 7850 | lambda=0.366 θ1=6.843 θ2=2.085 | acc(θ-block)=0.299
iter 7900 | lambda=0.479 θ1=6.838 θ2=2.094 | acc(θ-block)=0.299
iter 7950 | lambda=0.417 θ1=6.795 θ2=2.093 | acc(θ-block)=0.300
iter 8000 | lambda=0.437 θ1=6.837 θ2=2.081 | acc(θ-block)=0.300
iter 8050 | lambda=0.422 θ1=6.945 θ2=2.096 | acc(θ-block)=0.300
iter 8100 | lambda=0.424 θ1=6.814 θ2=2.091 | acc(θ-block)=0.299
iter 8150 | lambda=0.451 θ1=6.873 θ2=2.100 | acc(θ-block)=0.299
iter 8200 | lambda=0.439 θ1=6.843 θ2=2.085 | acc(θ-block)=0.299
iter 8250 | lambda=0.427 θ1=6.859 θ2=2.099 | acc(θ-block)=0.299
iter 8300 | lambda=0.430 θ1=6.893 θ2=2.094 | acc(θ-block)=0.298
iter 8350 | lambda=0.424 θ1=6.869 θ2=2.102 | acc(θ-block)=0.298
iter 8400 | lambda=0.390 θ1=6.865 θ2=2.085 | acc(θ-block)=0.298
iter 8450 | lambda=0.393 θ1=6.830 θ2=2.092 | acc(θ-block)=0.299
iter 8500 | lambda=0.430 θ1=6.809 θ2=2.077 | acc(θ-block)=0.299
iter 8550 | lambda=0.416 θ1=6.837 θ2=2.100 | acc(θ-block)=0.299
iter 8600 | lambda=0.402 θ1=6.826 θ2=2.085 | acc(θ-block)=0.300
iter 8650 | lambda=0.465 θ1=6.829 θ2=2.091 | acc(θ-block)=0.300
iter 8700 | lambda=0.415 θ1=6.816 θ2=2.096 | acc(θ-block)=0.300
iter 8750 | lambda=0.392 θ1=6.839 θ2=2.081 | acc(θ-block)=0.300
iter 8800 | lambda=0.433 θ1=6.839 θ2=2.099 | acc(θ-block)=0.299
iter 8850 | lambda=0.426 θ1=6.795 θ2=2.091 | acc(θ-block)=0.298
iter 8900 | lambda=0.430 θ1=6.912 θ2=2.083 | acc(θ-block)=0.298
iter 8950 | lambda=0.391 θ1=6.844 θ2=2.099 | acc(θ-block)=0.298
iter 9000 | lambda=0.433 θ1=6.830 θ2=2.102 | acc(θ-block)=0.298
iter 9050 | lambda=0.420 θ1=6.875 θ2=2.101 | acc(θ-block)=0.298
iter 9100 | lambda=0.401 θ1=6.876 θ2=2.093 | acc(θ-block)=0.297
iter 9150 | lambda=0.402 θ1=6.886 θ2=2.107 | acc(θ-block)=0.297
iter 9200 | lambda=0.428 θ1=6.846 θ2=2.092 | acc(θ-block)=0.296
iter 9250 | lambda=0.401 θ1=6.897 θ2=2.094 | acc(θ-block)=0.295
iter 9300 | lambda=0.413 θ1=6.791 θ2=2.093 | acc(θ-block)=0.295
iter 9350 | lambda=0.377 θ1=6.864 θ2=2.098 | acc(θ-block)=0.295
iter 9400 | lambda=0.418 θ1=6.830 θ2=2.100 | acc(θ-block)=0.296
iter 9450 | lambda=0.397 θ1=6.901 θ2=2.081 | acc(θ-block)=0.296
iter 9500 | lambda=0.363 θ1=6.867 θ2=2.092 | acc(θ-block)=0.295
iter 9550 | lambda=0.384 θ1=6.833 θ2=2.092 | acc(θ-block)=0.295
iter 9600 | lambda=0.405 θ1=6.873 θ2=2.087 | acc(θ-block)=0.295
iter 9650 | lambda=0.390 θ1=6.851 θ2=2.099 | acc(θ-block)=0.296
iter 9700 | lambda=0.423 θ1=6.871 θ2=2.081 | acc(θ-block)=0.296
iter 9750 | lambda=0.448 θ1=6.902 θ2=2.092 | acc(θ-block)=0.296
iter 9800 | lambda=0.417 θ1=6.867 θ2=2.097 | acc(θ-block)=0.296
iter 9850 | lambda=0.419 θ1=6.903 θ2=2.093 | acc(θ-block)=0.295
iter 9900 | lambda=0.435 θ1=6.832 θ2=2.102 | acc(θ-block)=0.295
iter 9950 | lambda=0.488 θ1=6.849 θ2=2.091 | acc(θ-block)=0.295
iter 10000 | lambda=0.423 θ1=6.845 θ2=2.100 | acc(θ-block)=0.295
save(fit1, file = "application_fit_hourly_final.RData")
df <- fit1$draws
g1 <- ggplot(df, aes(iter, theta1)) + geom_line() +
labs(title = expression(trace~theta[1]))
g2 <- ggplot(df, aes(iter, theta2)) + geom_line() +
labs(title = expression(trace~theta[2]))
g3 <- ggplot(df, aes(iter, lambda)) + geom_line() +
labs(title = expression(trace~lambda))
gd1 <- ggplot(df, aes(theta1)) + geom_density() +
labs(title = expression(density~theta[1]))
gd2 <- ggplot(df, aes(theta2)) + geom_density() +
labs(title = expression(density~theta[2]))
gd3 <- ggplot(df, aes(lambda)) + geom_density() +
labs(title = expression(density~lambda))
gridExtra::grid.arrange(g1,g2,g3, gd1,gd2,gd3, ncol = 3)
fit1$post_mean
lambda theta1 theta2
0.4186907 6.8560784 2.0937366
sd(fit1$draws$theta1)
[1] 0.03948329
sd(fit1$draws$theta2)
[1] 0.008061415
plot_extreme_on_map <- function(i, extreme_mat = extreme_L1_hourly,
map_sf = square_map, data_wide = data_hourly_sub) {
stopifnot(is.matrix(extreme_mat),
i >= 1, i <= nrow(extreme_mat))
# column order of pixels in the wide daily data used to build rain_mat/extreme_L1
pix_order <- as.integer(sub("P", "", colnames(data_wide)[-1]))
# values for the i-th extreme event
vals_i <- as.numeric(extreme_mat[i, ])
df_i <- data.frame(PIXEL = pix_order, val_i = vals_i)
# join to geometry
joined <- map_sf %>%
left_join(df_i, by = "PIXEL")
# quick stretch for nicer legend breaks
brks <- quantile(df_i$val_i, probs = seq(0, 1, length.out = 6), na.rm = TRUE)
brks_fixed <- sort(unique(brks))
# plot
mapview(
joined,
zcol = "val_i",
layer.name = paste0("Extreme L1 snapshot (row ", i, ")"),
at = brks_fixed
)
}
plot_extreme_on_map(1)
plot_extreme_on_map(5)
\(\chi^2(S_A, S_B)\) test
chi_q <- function(X, q = 0.98) {
n <- ncol(X)
m <- nrow(X)
# thresholds at site level
u <- apply(X, 2, quantile, probs = q, na.rm = TRUE)
# initialize matrix
chi_hat <- matrix(NA, n, n)
for (i in 1:(n-1)) {
for (j in (i+1):n) {
Ai <- X[, i] > u[i]
Aj <- X[, j] > u[j]
chi_hat[i, j] <- mean(Ai & Aj) / mean(Ai)
chi_hat[j, i] <- chi_hat[i, j]
}
}
diag(chi_hat) <- 1
chi_hat
}
dist_matrix <- as.matrix(dist(coord))
stopifnot(length(fit1$hard_labels) == nrow(extreme_L1_hourly))
labs <- fit1$hard_labels
X1 <- extreme_L1_hourly[labs == 1, , drop = FALSE]
X2 <- extreme_L1_hourly[labs == 2, , drop = FALSE]
# 4) compute χ̂_q for both comps and two q levels
chi_95_c1 <- chi_q(X1, q = 0.95)
chi_95_c2 <- chi_q(X2, q = 0.95)
chi_90_c1 <- chi_q(X1, q = 0.90)
chi_90_c2 <- chi_q(X2, q = 0.90)
upper_to_df <- function(chi_mat, dist_mat, q_label, comp_label) {
idx <- which(upper.tri(chi_mat), arr.ind = TRUE)
tibble(
dist = dist_mat[idx],
chi = chi_mat[idx],
q = q_label,
comp = comp_label
) %>% filter(is.finite(chi))
}
df_all <- bind_rows(
upper_to_df(chi_95_c1, dist_matrix, "q=0.95", "Comp1"),
upper_to_df(chi_95_c2, dist_matrix, "q=0.95", "Comp2"),
upper_to_df(chi_90_c1, dist_matrix, "q=0.90", "Comp1"),
upper_to_df(chi_90_c2, dist_matrix, "q=0.90", "Comp2"),
)
# 6) plot (mirrors your original style)
ggplot(df_all, aes(x = dist, y = chi, color = comp)) +
geom_point(alpha = 0.4, size = 0.6) +
facet_wrap(~ q, scales = "free_y") +
theme_bw() +
labs(
x = expression(h(s[A], s[B])),
y = expression(hat(chi)[q](s[A], s[B])),
title = "Empirical " %+% expression(hat(chi)[q]) %+% " vs Distance by MCMC Hard-Label Component"
) +
scale_color_manual(values = c("Comp1" = "blue", "Comp2" = "red"))
for (i in seq_len(ncol(fit1$z_draws))) {
plot(
fit1$z_draws[, i],
type = "l",
xlab = "iterations",
ylab = "z draw",
main = paste("sample", i) # <--- Modified line
)
}
plot(times_hourly, fit1$tau_bar, type = "l",
xlab = "Time", ylab = "Posterior mean of z",
main = "Posterior mean with season")